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Abstract 


Bovine leukemia virus (BLV) is the etiologic agent of enzootic bovine leucosis (EBL) for the bovine host. 
In this study to examine gene expression changes in the manifestation of the EBL malignancy, four pooled RNA 
samples (three RNAs in each sample) were applied for transcriptome sequencing using RNA-seq technique. 
Differential expression analysis was done to compare the infected bovine group with the healthy bovine group using 
DESeq2 package in R software. Furthermore, functional gene ontology (GO) term and KEGG pathway enrichment 
analysis were stablished using the DAVID online database to identify involved GO terms and pathways in the host 
response to BLV infection. Our results suggested that 371 up- and 72 downregulated genes were involved in EBL 
with statistically significant threshold log2foldchange (LFC) = | and false discovery rate (FDR) <0.05 that were 
enriched in 74 biological processes and 20 KEGG pathways. Most of identified genes were associated with cancer, 
especially B-cell malignancies. The glycolysis/glycogenesis metabolic process is activated in B cells that confers 
growth and survival advantages in tumor and dysregulated CXCL10, ILI7R, BTK, CDK4 and SYK genes known as 
valid biomarkers to increase the proliferation of malignant cell. The outcomes can provide a list of involved genes in 


the malignancy and help to screen candidate genes for cancer therapy in the future. 
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Introduction 

Bovine leukemia virus, a deltaretrovirus 
related to human T-cell leukemia virus types 1 and 2 
(HTLV-1 and HTLV-2), causes a lethal lymphoma 
or lymphosarcoma, enzootic bovine leucosis, which 
has been defined as a prolonged course that often 
involves persistent lymphocytosis (PL) and ends in 
B-cell lymphoma (Burny et al., 1988; Juliarena et al., 
2017). EBL happens in a small portion (< 5-10%) of 
infected adult bovine, while the majority of infected 
animals (~70%) remain asymptomatic carriers. The 
prevalence of BLV infection and inducing of EBL in 
carriers cause economic losses for the livestock 
industry in the worldwide (Polat et al., 2017) and this 
infection is as high as 25.4% in the northeast of Iran 
(Mousavi et al., 2014). 
The BLV virus carries two single-strand RNA 
(ssRNA), which can be reverse transcribed into 
DNA and integrated into the host genome as a 
proviral DNA (Moratorio et al., 2013). It encodes a 
trans-activator protein Tax, which plays a leading 
role in viral replication and activates dysregulation 
of cellular genes including cytokines, adhesion 
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molecules and growth factors at the initial steps of 
oncogenesis to affect viral spreading and disease 
progression (Arainga et al., 2012). Furthermore, Tax 
gene can transcriptionally change the expression 
level of cellular genes involved in the cell cycle, 
DNA repair and programmed cell death processes 
(Aida et al., 2013; Arainga et al., 2012; Kouznetsova 
et al., 2019) to provide the basis of cell 
transformation. Despite the activation of Tax 
protein, the expression of antisense transcripts AS1 
and AS2 from 3’LTR side of provirus have pivotal 
roles in the maintenance of malignancy and escaping 
from host immune responses, mainly at the latent 
stage of malignancy (Durkin et al., 2016). 

The viral infection and tax-responsive genes might 
not be unique supplier for the cell transformation 
while other additional events are required such as 
gene polymorphisms, chromosomal aberrations, 
genome instability and accumulation of mutations in 
cellular DNA (Klener et al., 2006) for progression of 
this multistep tumor in the bovine host. The 
existence of polymorphism in the enhancer region of 
TNF-a (Konnai et al., 2006) and mutation in 
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oncogenes such as P53 (Dequiedt et al., 1995) is 
confirmed to contribute in the progression of BLV- 
induced lymphoma. Therefore, the virus-host 
interaction manner apart from viral factors would 
have a principal role to induce malignancy. 

The objective of this study was to investigate the 
gene expression changes in the host-virus interaction 
manner using RNA-seq technique. It is known a little 
about the transcriptional changes associated with 
abnormal B-cell growth in EBL malignancy. We 
compared gene expression profiles of EBL and the 
healthy bovine groups to identify differentially 
expressed genes (DEGs). Subsequently, the 
identified DEGs were enriched to determine 
functional GO terms and KEGG pathways that 
participate in the B-cell transformation and EBL 
oncogenicity. 


Materials and Methods 


Sample collection and EBL diagnostics 

This study was approved and supervised by the 
ethics committee of Ferdowsi University of 
Mashhad, Mashhad, Iran (No. 42333). Six infected 
lymphoid tumors at the acute stage of the EBL 
disease, and six healthy lymph nodes without BLV 
viral infection collected from cows (3.5-4 years old), 
at the slaughterhouse of Mashhad, Iran, and 
immediately transferred to a liquid nitrogen tank. At 
first, the EBL disease approved by two veterinarians 
based on the existence of clinical symptoms. Then, a 
qualitative method was optimized using the Pol gene 
expression in Rotor-gene 6000 (Qiagen, Germany) 
machine. Reaction ingredients were composed of 10 
ul commercial Real-time master mix (Thermo 
Fisher, USA), 0.5 wl of primers (Fw, 5’- 
CCTCAATTCCCTTTAAACTA-3’; Rv, 5°- 
GTACCGGGAAGACTGGATTA-3’), 1 wuL of 
TaqMan probe (FAM- 
GAACGCCTCCAGGCCCTTCA-BHQ1) — (Rola- 
Luszczak et al., 2013), 2.5 pl of deionized water and 
5 ul of template DNA. DNAs were extracted with 
high quality by using the High Pure PCR template 
Kit (Qiagen, Germany), according to the 
manufacturer’s instruction. 


RNA extraction and sequencing 

Total RNAs were isolated using the RNeasy Mini 
Kit (Qiagen, Germany), according to the 
manufacturer’s instructions. The pooled RNAs per 
group (infected and uninfected) were prepared by 
mixing equal concentrations of three biological 
replicates (RNAs) to generate four pooled RNA 
samples. The single-strand DNA (ssDNA) templates 
were synthesized using oligo-dT primers and reverse 


http://jcmr.um.ac.ir 


15 


transcriptase, according to the manufacturer’s 
instructions (Thermo Fisher, USA). The purity, 
integrity and quantity of RNAs and cDNAs with 
amplification of the GAPDH gene (as a reference 
gene) were examined on 1.5% agarose gel and 
NanoDrop device. The type of poly-A capturing 
library was used for sequencing with Illumina HiSeq 
4000 platform to generate 2-100 bp paired-end 
reads, and 2-30 million raw paired-end read files, by 
the BGI Company. 


Differential expression analysis 

Raw reads were preprocessed using FastQC 
(Andrews, 2016) and Trimmomatic (Bolger et al., 
2014) software; then clean reads were acquired for 
subsequent analyses. The clean raw paired-end read 
profiles from each sample were mapped to the host- 
provirus reference genome individually, using 
STAR software (v2.5.4) (Dobin et al., 2013). Host- 
provirus reference genome was built by merging 
bovine reference genome UMD3.1 with proviral 
genome sequence BLV-YR2 (GenBank: KT 122858) 
before alignment. Sorting and indexing of STAR 
outputs were performed by SAMtools software 
(v1.9) (Li et al., 2009). Annotation of GTF files 
(v74) for bovine, was downloaded from Ensemble 
database, and custom BLV virus annotation was 
built by HTLV-1/BLV genomics lab in Medical 
Genomics/Unit of Animal Genomics (UAG), Liege, 
Belgium. The FeacureCounts package was used for 
reads quantification (Liao et al., 2013). Finally, 
differential expression analysis was performed 
based on statistically significant FDR < 0.05 in R 
software (v.3.5), using DESeq2 package (v.1.20) 
(Love et al., 2014). 


Gene ontology (GO) and pathway analysis 

Gene ontology analysis was performed on the list 
of significant (DEGs) based on biological process 
(BP), molecular function (MP), cellular component 
(CC) terms and the KEGG (The Kyoto Encyclopedia 
of Genes and Genomes) pathway, using DAVID 6.8 
web tools (https://david.ncifcrf.gov/) (Huang et al., 
2009). For identification of related GO terms to our 
samples, the GO profiles of significant DEGs were 
compared to the GO profile of bovine genome 
annotation as a reference set. The FDR < 0.05 and P- 
value < 0.01 were considered as a cutoff threshold 
for statistical significance in detection of GO terms 
and the KEGG pathways, respectively. The 
identified GO terms were classified using 
CateGOrizer “GO term classification counter” 
online tools based on GO slim method (Hu et al., 
2008). 
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Results 


EBL diagnostics 

The Enzootic bovine leucosis can be characterized 
by developing tumor masses, tumor masses which 
are composed of monomorphic lymphocytic cells, 
and variable clinical signs including cardiac 
lymphosarcoma, enlarged lymph nodes and discrete 
nodular masses or a diffuse tissue infiltrate in EBL 
samples that were shown in Figure 1. The 
amplification of Pol gene by real-time PCR 
technique confirmed the infection of EBL samples, 
while the healthy subjects were not seropositive for 
BLV infection (Figure 2). 


Figure 1. The clinical signs of EBL samples. a) Retro- 
ocular lymph node b) Cardiac lymphosarcoma, c) 
Enlarged mediastinal lymph node and d) Nodular 
masses. 
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Figure 2. RT-qPCR amplification curves of Pol gene 
in healthy and infected samples. 


RNA extraction and preprocessing of raw reads 
The appearance of 28S and 18S ribosomal RNA and 
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amplification of GAPDH gene as a reference on a 
1.5% agarose gel, had proven the quality of RNA 
and cDNA products (Figure 3). Pooled samples (2 
samples per group) were sequenced based on high 
throughput sequencing and by using of Illumia Hiseq 
4000 platform. Almost 2-30 million paired-end raw 
read profiles were obtained per sample from the 
company. All reads passed the preprocessing cutoff 
thresholds using FastQC software and only the 
adaptor contaminations were trimmed off by using 
Trimmomatic software. 


Figure 3. The integrity of RNA and complementary 
DNA amplification (CDNA). A) The RNA quality of 
samples based on 28S and 18S ribosomal RNA. B) 
cDNA quality with amplification of the GAPDH gene. 
Lane M shows 100 bp plus DNA marker. 


Differential expression analysis 

Read mapping was carried out using STAR 
aligner, a splice junction aware software. Most of the 
reads (60.7-79.5%) were mapped to the reference 
genome. The result of the DE analysis suggested 443 
unique genes that were differentially expressed in 
the infected group, compared to the healthy group. 
The identification of significant genes was based on 
the FDR <0.05. Around 371 genes were over- 
expressed with positive Log2—fold change more than 
2, and 72 genes were under-expressed with negative 
Log2 fold change less than 2, in infected samples 
versus control samples. The detailed information of 
the discovered genes is presented in Supplementary 
data file sheet1. 
Functional GO terms and pathway analysis 

The Gene Ontology enrichment analysis 

was accomplished for better understanding of the 
biological processes and involved pathways of 
identified DEGs in EBL malignancy. There were 
112 significant improved GO terms based on FDR < 
0.05 (30, 8 and 74 for CC, MF and BP categories, 
respectively) (see supplementary data file sheet 2). 
Because of redundancy and overlapping among 
identified GO term results, they were classified 
using GO slim method in CateGOrizer online tool to 
describe the basic terms and relationships between 
terms within the context of biology (Table 1). A term 
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Table 1. The classified Gene Ontology (GO) terms of DEGs, using CateGOrizer online tool. 


Parent GO ID | GO Definitions GO terms | Counts Fractions 
GO:0008150 | biological-process BP 74 36.63% 
GO:0008152 | metabolism BP 64 14.36% 
GO:0009058 | biosynthesis BP 29 11.39% 
GO:0006139 nucleobase, nucleoside, nucleotide and BP 3 1.49% 
nucleic acid metabolism 

GO:0016043 | cell organization and biogenesis BP 3 1.49% 
GO:0006259_ | DNA metabolism BP 3 0.99% 
GO:0006996 | organelle organization and biogenesis | BP 2 0.99% 
GO:0007049 | cell cycle BP 2 0.50% 
GO:0007165 | signal transduction BP 1 0.50% 
GO:0007154 | cell communication BP 1 14.36% 
GO:0005575 | cellular_component CC 29 41.43% 
GO:0005623 | cell CC 16 22.86% 
GO:0005622 | intracellular CC 13 18.57% 
GO:0005576 | extracellular region CC 5) 714% 
GO:0005634 | nucleus CC 4 5.71% 
GO:0005856 | cytoskeleton CC 1 1.43% 
GO:0005737 | cytoplasm CC 1 1.43% 
GO:0005730 | nucleolus CC 1 1.43% 
GO:0003674 | molecular_function MF 8 36.36% 
GO:0005488 | binding MF 8 36.36% 
GO:0003676 | nucleic acid binding MF 3 13.64% 
GO:0003723 | RNA binding MF 2 9.09% 
GO:0000166 | nucleotide binding MF 1 4.55% 


Only the basic (parent) GO terms are presented. BP, biological process; CC, cellular component; MF, molecular 
function; Counts are the number of identified GO terms in each class term; Fractions are the percentage of identified 


GO terms that grouped in each class term. 


Biological Process 


Cellular Component 


Molecular Function 


Figure 4. The classified basic Gene Ontology (GO) terms in pie graphs. 


may classify in a subset of different terms. The 
presented percentage in the fraction column of Table 
1, express the importance of the identified basic 
(parent) terms that consisted of child terms. The two 
top basic terms were metabolism and biosynthesis in 
BP category, cell and intracellular in CC category, 
binding and nucleic acid binding in MF category 
(Figure 4). Our result suggested 20 significant 
KEGG pathways based on statistically significant P- 
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value < 0.01 (Dahiru, 2008) (Table 2). The five top 
significant KEGG pathways were PI3K-Akt 
signaling pathway, pathway in cancer, biosynthesis 
of antibiotics, proteoglycans in cancer, cardiac 
muscle contraction and B cell receptor (BCR) 
signaling pathway, in EBL group. 
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Table 2. The enriched KEGG pathways of DEGs by using DAVID online tool 


KEGG ID Pathway Counts % Pop hits P-value 
bta04151 | PI3K-Akt signaling pathway 28 6.378 347 1.33E-05 
bta05200 | Pathways in cancer 30 6.833 398 2.18E-05 
bta01130 | Biosynthesis of antibiotics 20 4.555 206 2.74E-05 
bta05205 | Proteoglycans in cancer 19 4.328 203 7.49E-05 
bta04260 | Cardiac muscle contraction 11 2.505 80 2.04E-04 
bta04662 | B cell receptor signaling pathway 10 2.277 70 3.42E-04 
bta04916 | Melanogenesis 11 2.505 99 0.001 
bta04310 | Wnt signaling pathway 13 2.961 137 0.001 
bta05166 | HTLV-I infection 19 4.328 267 0.002 
bta04912 | GnRH signaling pathway 9 2.050 86 0.005 
bta04512 | ECM-receptor interaction 9 2.050 87 0.006 
bta04020 | Calcium signaling pathway 14 3.189 189 0.007 
bta05217 | Basal cell carcinoma 7 1.594 55 0.007 
bta05219 | Bladder cancer 6 1.366 40 0.008 
bta04915 | Estrogen signaling pathway 9 2.050 98 0.01 
bta00010 | Glycolysis / Gluconeogenesis 7 1.594 63 0.01 
bta04510 | Focal adhesion 14 3.189 208 0.01 
bta04660 | T cell receptor signaling pathway 9 2.05 105 0.01 
bta05414 | Dilated cardiomyopathy 8 1.82 86 0.01 
bta04611 | Platelet activation 10 2.27 127 0.01 


Only the significantly enriched KEGG pathways (P < 0.01) are presented. Counts are the number of unique DAVID 
gene symbol corresponding to our input gene list; Pop hits, the number of genes that belong to the specific gene 


category. 


Discussion 


The development and progression of 
leukemogenesis induced by BLV, depend on the 
level of viral factors such as Tax gene expression at 
the early stage of infection, while some recent 
studies have been reported that besides of infection, 
the immune deterioration, growth factor production 
and activation of cancer driver genes can induce cell 
transformation and malignancy to form acute phase 
of infection (Durkin et al., 2016; Gillet et al., 2007). 
Hence, the study of bovine’s transcriptome at the 
acute phase is helpful for a better understanding of 
involved genes and pathways in the manifestation of 
the EBL. 

In this study, we produced high throughput 
transcriptome data using RNA-seq technique from 4 
pooled RNA samples to examine gene expression 
changes in lymph node cells of two EBL samples in 
the counterpart of two sera-negative samples in 
response to BLV infection and progression to the 
lymphocytic stage. The expression of 443 unique 
DEGs was changed through BLV infection, 371 
were up- and 72 were downregulated. Of note, 
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previous studies were applied microarray technique 
(Brym and Kaminski, 2017; Everts et al., 2005; 
Klener et al., 2006) for investigation of gene 
expression changes, while we used RNA-seq 
technique for the first time. There are only a few 
DEGs in common between our result and the DEG 
lists of some microarray studies (e.g ADORA2B, 
CD19, APEX], HIFIA, CCT5, LMO2, HADH, 
DYNLLI, CIQBP, BCLIIA, SYK, PHGDH, 
WDFY4, NMI, CXXC5, SOX5, CYB5R3, RHNO1, 
YWHAG, LYN, Src, PEBP1, GNL3, CXXC5, ERH, 
AVEN, CMTM7, NACA, HSPA14, ENOPH1, LCN2, 
HMGB1, CXCL10, DYNLLI, CXHXorf57, ZBTB32, 
TUBB4B, and PKM) and a few DEGs with invert 
expression (e.g. HSPA4, ILI7R, CDK4, RNASE6, 
SDCBP, GAPDH, RGS2, MKI67, TOP2A, and 
CD79A) (Brym and Kaminski, 2017; Everts et al., 
2005; Klener et al., 2006). It is interesting that genes 
associated with several human malignancies, 
especially B-cell tumor, were represented among the 
differentially expressed genes. Some genes were 
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new and had no recognized role in tumor 
development. 
The upregulated CXCLIO was _ significantly 


associated with tumor cell migration and progression 
in diffuse large B-cell lymphoma (DLBCL) (Burny 
et al., 1988) and is a valid biomarker in progress of 
heart failure (Greenwood et al., 2018). In the EBL 
samples, the pathway of cardiac muscle contraction 
was significant, and heart failure might be a sign of 
clinical symptom of EBL malignancy. The 
upregulation of TOP2A and MKI67 in prostate 
cancer (Malhotra et al., 2011) and CDK4 in B cell 
lymphoma (Brym and Kaminski, 2017) are resulted 
to increase the proliferation of malignant cell and 
they might be good candidates for targeted cancer 
therapy. 

B cells are the main arm of humoral immunity, 
which is targeted by BLV, thus BCR and three 
membranes of B cell antigen receptors complex, 
CD79A, CD19, and CD70A, as most pivotal 
activators were upregulated in infected BLV cells 
(Brym and Kaminski, 2017). This strategy of BLV 
for B cell survival via BCR signaling (Burny et al., 
1988) is the homolog of the activation of TCR 
signaling in HTLV-1 infection (Ghezeldasht et al., 
2013) to induce cell cycle promoting proteins, 
provides growth factors in cell survival to inhibit 
apoptosis, and activates transforming pathways for 
malignancy. Activation of the BCR complex leads to 
the phosphorylation of proto-oncogenes tyrosine 
kinases family, including B7TK and SYK in 
hematopoietic cells (Gauld and Cambier, 2004). 
High expression of BTK and SYK genes increase the 
activation, proliferation, and survival of DLBCL 
(Gauld and Cambier, 2004; Tolar et al., 2005). 
Further, the identified differentially expressed 
genes were used in the gene enrichment analysis by 
using the online DAVID database to screen 
functional patterns associated with the expression 
differences observed. In general, our result 
suggested 74 BP, 30 CC and 8 MF. Because of the 
existence of more similarity and redundancy among 
obtained Go terms, we categorized GO terms in 
terms of basic GO classes using CateGorizer online 
tool based on GO slim method to reveal the main 
basic ‘parent? GO terms and the relationships 
between ‘child’ terms (Hu et al., 2008). Finally, 
whole expressed genes were localized in 7 cellular 
component term including cell, intracellular or 


extracellular region, nucleus, cytoskeleton, 
cytoplasm and nucleus; and were involved in 9 
biological processes including metabolism; 


biosynthesis; nucleobase, nucleoside, nucleotide and 
nucleic acid metabolism; cell organization and 
biogenesis; DNA metabolism; organelle 
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organization and biogenesis; cell cycle; signal 
transduction and cell communication which they 
mediated four molecular functions including 
binding; nucleic acid binding; RNA binding and 
nucleotide-binding. These enrichment results 
provided further supports to previously published 
results about the invasion of BLV infected cells 
(Brym and Kaminski, 2017; Everts et al., 2005; van 
den Heuvel et al., 2005). 

Furthermore, our study suggests the expression 
changes in genes from 20 distinct KEGG pathways, 
including PI3K-Akt signaling pathway; pathway in 
cancer; biosynthesis of antibiotics; proteoglycans in 
cancer; cardiac muscle contraction; B cell receptor 
signaling pathway; melanogenesis; Wnt signaling 


pathway; HTLV-1 infection; GnRH _ signaling 
pathway; ECM-receptor interaction; calcium 
signaling pathway; Basal cell carcinoma; Bladder 
cancer; estrogen signaling pathway; 


glycolysis/gluconeogenesis; focal adhesion; T-cell 
receptor signaling pathway; dilated cardiomyopathy 
and platelet activation. The dysregulation and 
activation of mentioned KEGG pathways, leading to 
the disruption of B-cell survival and growth control. 
Interestingly, these activated processes were 
reported previously with HTLV-1, a T-cell-tropic 
oncoretrovirus, as well as non-virus-associated 
human B-cell malignancies (Liu et al., 2018; 
Mozhgani et al., 2018). 

The activation of cellular adhesion molecules 
(CAMs) and focal adhesion process play central 
roles in transferring and localization of leukemic 
cells in different tissues that finally affect cancer cell 
growth and metastasis (Brym and Kaminski, 2017; 
Burny et al., 1988). More studies are needed to 
explain the impact of the adhesion process in BLV 
infiltration to host cells. The 
glycolysis/gluconeogenesis pathway happens inside 
the cells to produce adenosine triphosphate (ATP), 
as a cell energy, to drive many processes in cells. 
Recent studies have reported the high activation of 
glycolysis/gluconeogenesis, PI3K-AKT signaling 
pathway, B cell and T cell receptor signaling 
pathway in cancer cells, including colorectal cancer, 
adult T cell leukemia (ATLL), and diffuse large B 
cell lymphoma (DLBCL) compared to normal cells 
(Caro et al., 2012; Liu et al., 2018; Mozhgani et al., 
2018). Our results indicated that the glycolysis 
metabolic process was activated in B cells that 
confers growth and survival advantages in the tumor. 


In conclusion, a wide range of effective, up- and 
downregulated genes involved in EBL pathogenicity 
and malignancy were introduced by the present 
study. More identified processes including ECM- 
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receptor interaction, focal adhesion, BCR signaling, 
TCR signaling, glycolysis/gluconeogenesis, PI3K- 
Akt pathway and the dysregulated genes including 
CXCLIO, BTK, SYK, ILI7R and CDK4 were 
associated with the disruption of natural DNA 
replication and cell proliferation which induce the 
repression of DNA damage responses and apoptosis 
in B-cells to make malignancy. Moreover, the over- 
activation of glycolysis/gluconeogenesis pathway 
along with immune cells is required to support the 
metabolic reprogramming of B cells. As regards 
EBL malignancy is a life-threatening disease for 
farm animals, more researches on transcriptome 
level are necessary to accelerate our findings of the 
EBL diseases. 
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